Normalization and Batch Effect Correction: Multiome
28 May 2021
1 Objective
The objective of this notebook is to perform pre-processing and dimensional reduction on both assays independently, using standard approaches for RNA and ATAC-seq data. Then, we will follow the “Joint RNA and ATAC analysis: 10x multiomic” vignette from Signac to obtain a joint visualization using both modalities.
2 Pre-processing
2.1 Load packages
library(Signac)
library(Seurat)
library(GenomicRanges)
library(future)
library(SeuratWrappers)
library(harmony)
library(EnsDb.Hsapiens.v86)
library(stringr)
library(dplyr)
library(ggplot2)
set.seed(173)2.2 Parameters
path_to_obj <- here::here("multiome/results/R_objects/4.tonsil_multiome_filtered_combined_with_metadata.rds")
path_to_save <- here::here("multiome/results/R_objects/5.tonsil_multiome_integrated_using_wnn.rds")2.3 Load Multiome filtered data
tonsil_filtered <- readRDS(path_to_obj)
tonsil_filtered## An object of class Seurat
## 188447 features across 47150 samples within 2 assays
## Active assay: peaks (151846 features, 0 variable features)
## 1 other assay present: RNA3 Without Harmony integration
3.1 scATAC
3.1.1 Normalization and linear dimensional reduction
DefaultAssay(tonsil_filtered) <- "peaks"
tonsil_filtered <- RunTFIDF(tonsil_filtered)
tonsil_filtered <- FindTopFeatures(tonsil_filtered, min.cutoff = "q0")
tonsil_filtered <- RunSVD(tonsil_filtered)DepthCor(tonsil_filtered)
3.1.2 UMAP representation
tonsil_filtered <- RunUMAP(
tonsil_filtered,
dims = 2:40,
reduction = "lsi",
reduction.name = "umap.atac",
reduction.key = "atacUMAP_"
)
DimPlot(
tonsil_filtered,
reduction = "umap.atac",
group.by = "library_name",
pt.size = 0.1
) + NoLegend()
3.2 scRNA
3.2.1 Normalization and linear dimensional reduction
DefaultAssay(tonsil_filtered) <- "RNA"
tonsil_filtered <- NormalizeData(
tonsil_filtered,
normalization.method = "LogNormalize",
scale.factor = 1e4
)
tonsil_filtered <- tonsil_filtered %>%
FindVariableFeatures(nfeatures = 3000) %>%
ScaleData() %>%
RunPCA() 3.2.2 UMAP representation
tonsil_filtered <- RunUMAP(
tonsil_filtered,
dims = 1:30,
reduction = "pca",
reduction.name = "umap.rna",
reduction.key = "rnaUMAP_"
)
DimPlot(
tonsil_filtered,
reduction = "umap.rna",
group.by = "library_name",
pt.size = 0.1
) + NoLegend()
3.3 Joint
#tonsil_filtered <- FindMultiModalNeighbors(
# tonsil_filtered,
# reduction.list = list("pca", "lsi"),
# dims.list = list(1:30, 2:40)
#)
# tonsil_filtered <- RunUMAP(
# tonsil_filtered,
# nn.name = "weighted.nn",
# reduction.name = "wnn.umap",
# reduction.key = "wnnUMAP_"
# )
# DimPlot(
# tonsil_filtered,
# reduction = "wnn.umap",
# group.by = "library_name",
# pt.size = 0.1
# ) + NoLegend()4 With Harmony integration
4.1 scATAC
DefaultAssay(tonsil_filtered) <- "peaks"
tonsil_filtered <- RunHarmony(
object = tonsil_filtered,
reduction = "lsi",
dims = 2:40,
group.by.vars = "gem_id",
assay.use = "peaks",
project.dim = FALSE,
reduction.save = "harmony_peaks"
)4.1.1 UMAP representation
tonsil_filtered <- RunUMAP(
tonsil_filtered,
dims = 2:40,
reduction = "harmony_peaks",
reduction.name = "umap.atac",
reduction.key = "atacUMAP_"
)
DimPlot(
tonsil_filtered,
reduction = "umap.atac",
group.by = "library_name",
pt.size = 0.1
) + NoLegend()
4.2 scRNA
DefaultAssay(tonsil_filtered) <- "RNA"
tonsil_filtered <- RunHarmony(
object = tonsil_filtered,
reduction = "pca",
dims = 1:30,
group.by.vars = "gem_id",
assay.use = "RNA",
project.dim = FALSE,
reduction.save = "harmony_RNA"
)tonsil_filtered <- RunUMAP(
tonsil_filtered,
dims = 1:30,
reduction = "harmony_RNA",
reduction.name = "umap.rna",
reduction.key = "rnaUMAP_"
)
DimPlot(
tonsil_filtered,
reduction = "umap.rna",
group.by = "library_name",
pt.size = 0.1
) + NoLegend()
4.3 Joint
# tonsil_filtered <- FindMultiModalNeighbors(
# tonsil_filtered,
# reduction.list = list("harmony_peaks", "harmony_RNA"),
# dims.list = list(2:40, 1:30)
# )
# tonsil_filtered <- RunUMAP(
# tonsil_filtered,
# nn.name = "weighted.nn",
# reduction.name = "wnn.umap",
# reduction.key = "wnnUMAP_"
# )
# DimPlot(
# tonsil_filtered,
# reduction = "wnn.umap",
# group.by = "library_name",
# pt.size = 0.1
# ) + NoLegend()5 Save
We will save the resulting object and use it in the following notebook to exclude doublets:
saveRDS(tonsil_filtered, path_to_save)6 Session Information
sessionInfo()## R version 4.0.3 (2020-10-10)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: macOS Big Sur 10.16
##
## Matrix products: default
## BLAS/LAPACK: /Users/pauli/opt/anaconda3/envs/Tonsil_atlas/lib/libopenblasp-r0.3.10.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggplot2_3.3.2 dplyr_1.0.2 stringr_1.4.0 EnsDb.Hsapiens.v86_2.99.0 ensembldb_2.12.1 AnnotationFilter_1.12.0 GenomicFeatures_1.40.1 AnnotationDbi_1.50.3 Biobase_2.48.0 harmony_1.0 Rcpp_1.0.5 SeuratWrappers_0.3.0 future_1.20.1 GenomicRanges_1.40.0 GenomeInfoDb_1.24.0 IRanges_2.22.1 S4Vectors_0.26.0 BiocGenerics_0.34.0 Seurat_3.9.9.9010 Signac_1.1.0.9000 BiocStyle_2.16.1
##
## loaded via a namespace (and not attached):
## [1] reticulate_1.18 tidyselect_1.1.0 RSQLite_2.2.1 htmlwidgets_1.5.2 grid_4.0.3 BiocParallel_1.22.0 Rtsne_0.15 munsell_0.5.0 codetools_0.2-17 ica_1.0-2 miniUI_0.1.1.1 withr_2.3.0 colorspace_2.0-0 OrganismDbi_1.30.0 knitr_1.30 rstudioapi_0.12 ROCR_1.0-11 tensor_1.5 listenv_0.8.0 labeling_0.4.2 GenomeInfoDbData_1.2.3 polyclip_1.10-0 farver_2.0.3 bit64_4.0.5 rprojroot_2.0.2 parallelly_1.21.0 vctrs_0.3.4 generics_0.1.0 xfun_0.18 biovizBase_1.36.0 BiocFileCache_1.12.1 lsa_0.73.2 ggseqlogo_0.1 R6_2.5.0 rsvd_1.0.3 bitops_1.0-6 spatstat.utils_1.17-0 reshape_0.8.8 DelayedArray_0.14.0 assertthat_0.2.1 promises_1.1.1 scales_1.1.1
## [43] nnet_7.3-14 gtable_0.3.0 globals_0.13.1 goftest_1.2-2 ggbio_1.36.0 rlang_0.4.8 RcppRoll_0.3.0 splines_4.0.3 rtracklayer_1.48.0 lazyeval_0.2.2 dichromat_2.0-0 checkmate_2.0.0 BiocManager_1.30.10 yaml_2.2.1 reshape2_1.4.4 abind_1.4-5 backports_1.2.0 httpuv_1.5.4 Hmisc_4.4-1 RBGL_1.64.0 tools_4.0.3 bookdown_0.21 ellipsis_0.3.1 RColorBrewer_1.1-2 ggridges_0.5.2 plyr_1.8.6 base64enc_0.1-3 progress_1.2.2 zlibbioc_1.34.0 purrr_0.3.4 RCurl_1.98-1.2 prettyunits_1.1.1 rpart_4.1-15 openssl_1.4.3 deldir_0.2-3 pbapply_1.4-3 cowplot_1.1.0 zoo_1.8-8 SummarizedExperiment_1.18.1 ggrepel_0.8.2 cluster_2.1.0 here_1.0.1
## [85] magrittr_1.5 RSpectra_0.16-0 data.table_1.13.2 lmtest_0.9-38 RANN_2.6.1 SnowballC_0.7.0 ProtGenerics_1.20.0 fitdistrplus_1.1-1 matrixStats_0.57.0 hms_0.5.3 patchwork_1.1.0 mime_0.9 evaluate_0.14 xtable_1.8-4 XML_3.99-0.3 jpeg_0.1-8.1 gridExtra_2.3 compiler_4.0.3 biomaRt_2.44.4 tibble_3.0.4 KernSmooth_2.23-17 crayon_1.3.4 htmltools_0.5.0 mgcv_1.8-33 later_1.1.0.1 Formula_1.2-4 tidyr_1.1.2 DBI_1.1.0 tweenr_1.0.1 dbplyr_1.4.4 MASS_7.3-53 rappdirs_0.3.1 Matrix_1.2-18 igraph_1.2.6 pkgconfig_2.0.3 GenomicAlignments_1.24.0 foreign_0.8-80 plotly_4.9.2.1 xml2_1.3.2 XVector_0.28.0 VariantAnnotation_1.34.0 digest_0.6.27
## [127] sctransform_0.3.1 RcppAnnoy_0.0.16 graph_1.66.0 spatstat.data_1.4-3 Biostrings_2.56.0 rmarkdown_2.5 leiden_0.3.5 fastmatch_1.1-0 htmlTable_2.1.0 uwot_0.1.8.9001 curl_4.3 shiny_1.5.0 Rsamtools_2.4.0 lifecycle_0.2.0 nlme_3.1-150 jsonlite_1.7.1 viridisLite_0.3.0 askpass_1.1 BSgenome_1.56.0 pillar_1.4.6 lattice_0.20-41 GGally_2.0.0 fastmap_1.0.1 httr_1.4.2 survival_3.2-7 remotes_2.2.0 glue_1.4.2 spatstat_1.64-1 png_0.1-7 bit_4.0.4 ggforce_0.3.2 stringi_1.5.3 blob_1.2.1 latticeExtra_0.6-29 memoise_1.1.0 irlba_2.3.3 future.apply_1.6.0